# Back-ward covariable selection based on p-value
#
# @fittedModel: object reterned from function 'agModeling'
# @alpha: significant level
# Author: guochun
###############################################################################


mutliBackwardSelecter=function(fittedModel,alpha=0.05,initalnu=0.5){
	start=TRUE
	dels=list(NA)
	dels0=NA
	LastAIC=AIC(attr(fittedModel,"re.ppm"))
	currentAIC=0
	while(start | !is.na((dels[[1]])[1])){
		start=FALSE
		if(currentAIC<=LastAIC){
			model0=fittedModel
			dels0=dels
			if(currentAIC!=0)
				LastAIC=currentAIC
			dels=multiVariableSelection(fittedModel,alpha)
			#if there is no non-significant varable existed
			if(is.na(dels[[1]]))
				dels0=dels
		}else{
			if(is.na(dels0))
				dels0=dels
			dels=list(NA)
			fittedModel=model0
		}
		if(!is.na(dels[[1]])){
			tempAIC=numeric()
			tempmodel=list()
			for (i in 1:length(dels)){
				del=dels[[i]]
				tempi=constructModel(fittedModel@population, fittedModel@covr[-del],
						select=fittedModel@select, pvalue=fittedModel@pvalue, models=list(attr(fittedModel,"models")),
						nu=initalnu)
				tempmodel[[i]]=agModeling(tempi)
				tempAIC[i]=AIC(attr(tempmodel[[i]],"re.ppm"))
			}
			imin=which(tempAIC==min(tempAIC))[1]
			fittedModel=tempmodel[[imin]]
			currentAIC=tempAIC[imin]
			ap=attr(dels,"pvalues")
			dels=dels[imin]
			attr(dels,"pvalues")=ap
		}
		
	}
	del=dels0
	
	#calculate global pvalue
	if(length(fittedModel@covr)!=0 ){
		#print(del)
		pForH=prod(Sidak(attr(del,"pvalues")))
	}else{
		pForH=1
	}
	attr(fittedModel,"global_pForH")=pForH
	trend=predict(attr(fittedModel,"re.ppm"),type="trend")
	attr(fittedModel,"trend")=trend
	return(fittedModel)
}

AICbackwardSelecter=function(fittedModel,alpha=0.05,initalnu=0.5){
	start=TRUE
	del=NA
	del0=NA
	LastAIC=AIC(attr(fittedModel,"re.ppm"))
	currentAIC=0
	while(start | !is.na(del)){
		start=FALSE
		if(currentAIC<=LastAIC){
			model0=fittedModel
			del0=del
			if(currentAIC!=0)
			    LastAIC=currentAIC
			del=variableSelection(fittedModel,alpha)
		}else{
			del0=del
			del=NA
			fittedModel=model0
		}
		if(!is.na(del)){
			fittedModel=constructModel(fittedModel@population, fittedModel@covr[-del],
					select=fittedModel@select, pvalue=fittedModel@pvalue, models=list(attr(fittedModel,"models")),
					nu=initalnu)
			fittedModel=agModeling(fittedModel)
			currentAIC=AIC(attr(fittedModel,"re.ppm"))
		}
		
	}
	del=del0
	
	#calculate global pvalue
	if(length(fittedModel@covr)!=0 ){
		pForH=prod(Sidak(attr(del,"pvalues")))
	}else{
		pForH=1
	}
	attr(fittedModel,"global_pForH")=pForH
	trend=predict(attr(fittedModel,"re.ppm"),type="trend")
	attr(fittedModel,"trend")=trend
	return(fittedModel)
}



backwardSelecter=function(fittedModel,alpha=0.05,initalnu=0.5){
	start=TRUE
	del=NA
	while(start | !is.na(del)){
		start=FALSE
		del=variableSelection(fittedModel,alpha)
		if(!is.na(del)){
			fittedModel=constructModel(fittedModel@population, fittedModel@covr[-del],
					select=fittedModel@select, pvalue=fittedModel@pvalue, models=list(attr(fittedModel,"models")),
					nu=initalnu)
			fittedModel=agModeling(fittedModel)
		}
	}
	
	#calculate global pvalue
	if(length(fittedModel@covr)!=0 ){
		pForH=prod(Sidak(attr(del,"pvalues")))
	}else{
		pForH=1
	}
	attr(fittedModel,"global_pForH")=pForH
	trend=predict(attr(fittedModel,"re.ppm"),type="trend")
	attr(fittedModel,"trend")=trend
	return(fittedModel)
}



Sidak <- function(vecP)
#
# This function corrects a vector of probabilities for multiple testing
# using the Bonferroni (1935) and Sidak (1967) corrections.

{
	k = length(vecP)
	
	vecPB = 0
	vecPS = 0
	
	for(i in 1:k) {
		bonf = vecP[i]*k
		if(bonf > 1) bonf=1
		vecPB = c(vecPB, bonf)
		vecPS = c(vecPS, (1-(1-vecP[i])^k))
	}
#
	return(SidakP=vecPS[-1])
}


